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Abstract 

A new numerical method is developed for the solution of the two-dimensional, steady 
Navier-Stokes equations. The method that is presented differs in significant ways from 
the established numerical methods for solving the Navier-Stokes equations. The major 
differences are the following: First, the focus of the present method is on satisfying flux 
conservation in an integral formulation, rather than on simulating conservation laws in their 
differential form. Second, the present approach provides a unified treatment of the discrete 
dependent variables and their derivatives. All are treated as unknowns together to be 
solved for through simulating local and global flux conservation. Third, fluxes are balanced 
at cell interfaces without the use of interpolation or flux limiters. Fourth, flux conservation 
is achieved through the use of discrete regions known as conservation elements and solution 
elements. These elements are not the same as the standard control volumes used in the 
finite-volume method. Fifth, the discrete approximation obtained on each solution element 
is a functional solution of both the integral and differential form of the Navier-Stokes 
equations. Finally, the method that is presented is a highly localized approach in which 
the coupling to nearby cells is only in one direction for each spatial coordinate, and involves 
only the immediately adjacent cells. 

A general third-order formulation for the steady, compressible Navier-Stokes equa- 
tions is presented, and then a Newton’s method scheme is developed for the solution of 
incompressible, laminar channel flow. It is shown that the Jacobian matrix is nearly block 
diagonal if the nonlinear system of discrete equations is arranged appropriately and a 
proper pivoting strategy is used. Numerical results are presented for Reynolds numbers of 
100, 1000, and 2000. Finally, it is shown that the present scheme can resolve the developing 
channel flow boundary layer using as few as six to ten cells per channel width, depending 
on the Reynolds number. 




A New Flux Conserving Newton’s Method Scheme for the 
Two-Dimensional, Steady Navier-Stokes Equations 

I. Introduction 

This paper is concerned with the development of a new numerical approach, called 
the space-time solution element method,* for solving the Navier-Stokes equations. The 
present work builds on the ideas recently presented by Chang and To 1,2 for the numerical 
solution of conservation laws, and is part of a larger effort to build from the ground up a 
new family of flux conserving numerical schemes for solving the unsteady Navier-Stokes 
equations. 

The numerical framework that we are proposing differs in significant ways from the 
well established traditional methods for solving the Navier-Stokes equations. The major 
differences are the following: First, our approach is based on satisfying flux conservation 
in a space-time integral formulation. The focus is not on simulating partial differential 
equations in their differential form, but rather in satisfying space-time flux conservation 
of physical quantities. Second, the method we propose is a unified approach, in that it 
provides a unified treatment of space and time, and a unified treatment of the discrete 
dependent variables and their derivatives. All are treated as unknowns together to be 
solved for through simulating local and global flux conservation. Third, fluxes are balanced 
at cell interfaces without interpolation or the use of flux limiters. This differs fundamentally 
from a typical finite- volume method. 3 Fourth, local and global flux conservation is achieved 
through the use of discrete regions in space-time known as conservation elements and 
solution elements. These elements are not the same as the standard control volumes 
used in the finite-volume method. Fifth, the numerical approximation obtained on each 
solution element is a functional solution which satisfies the Navier-Stokes equations in 
both integral and differential form to some specified order. Finally, our method is a highly 
localized approach. The coupling to nearby cells is only in one direction for each space-time 
coordinate (as in a first-order method) and involves only the immediately adjacent cells. 
This feature holds irrespective of the order of accuracy of a particular numerical scheme 
based on our general approach. 

To better illustrate these distinguishing features, let 

V • h = 0 (1.1) 

be one of the governing equations of fluid mechanics, where the time derivative is part of 
the divergence operator, and h is a space-time flux current density vector. Using Gauss’ 
divergence theorem, equation (1.1) may be recast in an equivalent space-time integral form 

(j> h ■ ds = 0. (1.2) 

Js(V) 

*An abbreviation for The Method of Space-Time Conservation Element and Solution 
Element. 
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5(F) is the surface inclosing an arbitrary volume F in space-time, and ds is equal to dan 
where n is the outward unit normal to 5(F) and da is the area of a surface element of 
5(F). 

Equation (1.2) expresses a physical conservation law, namely, that the total flux out of 
5(F) of h is equal to zero. In general, equation (1.2) implies that there is a physical quantity 
which is conserved in the space-time volume F. The authors would like to emphasize that, 
unlike equation (1.1), equation (1.2) remains valid even when the components of h become 
discontinuous (Indeed, equation (1.1) is derived from equation (1.2) assuming that h is 
smooth). This feature, coupled with the fact that equation (1.2) expresses a physical 
conservation law which is valid for any space-time volume F, suggests that, in general, 
equation (1.2) is a better starting point for numerical calculations than is equation (1.1). 
The question then becomes how to best numerically simulate equation (1.2). 

The most general and natural way to represent the exact solution h by an approximate 
solution h in some discrete region F in space-time is through a Taylor series expansion. 

The discrete analogue to equation (1.2) would then be 

<£ h-d& = 0. (1.3) 

Js(V) ~ 

Note that to faithfully represent equation (1.2) in a discrete way, equation (1.3) must be 
satisfied on every computational cell and every union of computational cells in a mesh 
which has been used to discretize a flow field. Equation (1.3) is thus the starting point for 
any numerical algorithm based on the space-time solution element method. 

In this paper we are concerned with applying the ideas outlined above to the steady 
Navier-Stokes equations. Our ultimate aim, as expressed above, is to develop a family of 
new flux conserving numerical schemes for the unsteady Navier-Stokes equations. To lay 
the foundation for this effort, the present work will concentrate on the steady equations. 

The specific application that we address in this paper is incompressible, laminar chan- 
nel flow. A flux conserving Newton’s method scheme is developed based on the ideas 
expressed above and applied to the channel flow problem. The numerical method that 
we present differs in fundamental ways from any other Newton’s method scheme that the 
authors are aware of. Previous work for the Navier-Stokes equations has typically involved 
a finite-difference formulation. References 4 - 10 are a partial listing of work by other 
authors. 

We begin the development of our scheme in the next section by formulating the steady, 
compressible Navier-Stokes equations in the form of conservation laws. Following that, 
we derive the discrete flux conservation equations for a mesh with coinciding rectangular 
conservation elements and solution elements. We then use the discrete equations to develop 
an implicit scheme for solving incompressible, laminar channel flow. Finally, we present 
and discuss numerical results. 
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II. Conservation Laws for the Navier-Stokes Equations 


We consider the two-dimensional, steady, compressible Navier-Stokes equations in 
dimensionless form. ^Ve assume that the ratio of specific heats 'y, the viscocity //, and the 
coefficient of thermal conductivity k are all constant, and that the fluid is a peifect gas. 
Let Rcl, — P 00 ^ 00 ^ denote the Reynolds number, let Pr = denote the Prandtl number 
(where C p is the specific heat at constant pressure), and let denote the 

free-stream Mach number (where R is a gas constant). The parameters L, Uoot P<x>i and 
Too refer to some reference length, velocity, density, and temperature, respectively. 

Let x and y denote the horizontal and vertical coordinates, respectively, of a two- 
dimensional Euclidean space Denoting the fluid density by p, the horizontal velocity 
component by u, the vertical velocity component by u, the static pressure by p, and the 
temperature by T, the governing continuity, momentum, and energy equations may then 
be written in Cartesian coordinates as 11 


d(pu) + d(pv) _ Q 


dx 


dy 


£(,** + p - T xx ) + -J -( pUV - T xy ) = 0 

OX dy 

d d 

— — T xy ) + ‘q^(p v + P — T yy ) = ^ 


dx 


{[,< 


e + 


u 2 + v 2 


) + p 


'll U T x x V T x y 


} 


Off, U 2 + v 2 ] 1 

+ \p( e H 2 ) + P V — U T X y V Tyy + ? j, j ~ 


dy 


( 2 . 1 ) 

(2.2) 

(2.3) 


(2.4) 


where 


T xy — 


' yy 


2 (2 du _ dv, 

3 Re L dx dy 
1 du dv 

Re L dy dx 

2 (2 dv _ du 

3 Re,}" dy dx’ 


e = 


M 2 ^ 7(7-1) 


(2.5) 


(2.6) 

(2.7) 

( 2 . 8 ) 
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(2.9) 


1 dT 

(7 - 1 )Ml 0 Re L Pr Hi 

1 dT 

(7 - 1 )Ml 0 Re L Pr dy 


( 2 . 10 ) 


By applying Gauss’ divergence theorem to equations (2.1) - (2.4), they may be written 
in integral form as 

<fi h M -ds= 0 ( 2 . 11 ) 

Js(V) 

Y hxM ' ds = 0 (2.12) 

Js(V) 

Y hyM ' ds = 0 (2.13) 

JS(V) 

<p h E -ds= 0 (2.14) 

JS(V) 

where S(V) is the boundary of an arbitrary region V in E 2 . The flux current den- 
sity vectors, h M , h XM , h YM , and h E , corresponding to the continuity, x-momentum, y- 
momentum, and and energy equations, respectively, are defined by 

h M d = (pu,pv) (2.15) 


7* def / 9 

h X M = {pu + p - T xx , pUV - T xy ) (2.16) 


def 


hyM “ (pUV — T XJ/ , pV + P — Tyy) 


(2.17) 


- (h- 


u 2 + V 2 


e-\ ) + p u - UT XX - vr xy + q x , 


’ u 2 + V 2 ] 1 

P( e + 2 / + P^v - UT xy - V Tyy + q y > 


(2.13) 


Equations (2.11) - (2.14) thus represent physical conservation laws for the conservation 
of mass, momentum, and energy in an arbritrary region V of E 2 . 
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III. Derivation of the Discrete Flux Conservation Equations 


Let E 2 be discretized by a mesh with nonoverlapping rectangular regions. We assume 
constant spacing Ax in the x direction and constant spacing Ay in the y direction. (See 
Figure 1.) Each of the rectangular regions in the mesh will be referred to as both a 
conservation element and a solution element of E 2 . A conservation element is a discrete 
region in E 2 over which the flux conservation constraints (2.11) - (2.14) will be imposed. 
A solution element is a discrete region in E 2 in which a local Taylor series expansion will 
be employed to represent the physical solution. In general, they need not refer to the 
same discrete region. (See [2] for an example where the conservation elements and solution 
elements do not refer to the same discrete region.) A conservation element will be denoted 
by CE(i,j ) and a solution element by SE(i,j). The boundary of a conservation element 
will be denoted by S(CE(i,j)). We will denote the cell centers of the conservation and 
solution elements by ( Xi,y } ), where the subscript i is an index for the x coordinates, and 
the subscript j is an index for the y coordinates. 

We then assume that the fluid density, u and v velocity, pressure, and temperature 
can each be represented locally on a solution element by a two-dimensional Taylor series 
expansion about the cell center (xj,yj) as follows: 

p(x,y;i,j) = pi,o( x ~ x i) 2 + PiA x ~ x i)(v ~ yj) + Po.-iiy ~ Vj) 2 

+ p lt 0 (x - x,) + p 0 ,i(y - yj) + po.o (3.1) 

u(x,y;i,j) d = u 2>0 (x - x,) 2 + u ltl (x - x,)(y - yj) + u 0 , 2 (y - yj) 2 

+ u uo (x - Xi) + u 0fl (y - y } ) + u 0 0 (3.2) 


u(x,y;i,j) = v 2 ' 0 (x - Xi ) 2 + v ltl (x - Xi)(y - yj) + v 0i2 (y - yj) 2 

+ v lj0 (x - Xi) + v 0 il (y -yj) + v 0> 0 (3.3) 

p{x, y; i,j) d = p 2 , 0 ( x - x t) 2 +PiA x - x i)(y -yj) + p<>Ay -yj ) 2 

+ Pj 0 (x - Xi) + Po,,(y - y } ) + Po.o (3.4) 

T(x,y;i,j) d = T 2 , 0 (x - x,) 2 + T t>1 (x - Xi)(y - y 3 ) + T 0 , 2 (y - y ; ) 2 

+ T, 0 (x — x,) + T 0 ,i(y~yj) + T o.o (3-5) 

For clarity the i,j subscripts have been omitted from the coefficients in the Taylor series 
expansions. These coefficients are the unknowns to be solved for. 
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With this notation, the Taylor series coefficients are related to the derivatives of the 
discrete dependent variables at the cell center as follows: 

u 0 , o = u 
u l0 = du/dx 
u 0 , = du/dy 

u 2 ,o = ^ d 2 u/dx 2 

U 0 , 2 = i d 2 u/dy 2 
Uj , = d 2 u/dxdy 

and similarly for p,v,p, and T. 

The discrete analogue to equations (2.11) - (2.14) over an arbitrary conservation ele- 
ment CE(i,j ) in E 2 is then given by 


® h M • ds = 0 

JS{CE(i,j)) ~ 

(3.6) 

y h XM • ds — 0 

JS(CE(i,j)) 

(3.7) 

O 

II 

t <*> 

* 

sr 

0 

(3.8) 

<p h E • ds = 0 

JS(CE(i,j)) ~ 

(3.9) 

where 


h M d = (pu,pv) 

(310) 

d e f / 2 , \ 

= (PU + P ~ T xx , pUV — T xy ) 

(3.11) 
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(3.12) 


hvM — Txyi P Z P Zyy) 


h E d = { [/»(c + (u 2 + v 2 )/2^ + p 


w w T ri ZZ x y Q* 


P (e + (U 2 + U 2 )/2) + P U - U T X y - VTyy + (} y J (3.13) 


and 


Zxz = -^—{‘Zdu/dx - dy/dy ) (3-14) 

OilCj; 

T xy = 4-{du/dy + du/ds) (3.15) 

Iic L 

Zyy = rj— (2dv/dy - du/dx) (3.16) 

0/16^ 


e 


i 

M£,7(7-l) 


T 


(3.17) 


9ar 


?» 


l 

( 7 - 


dT/dx 


1 

( 7-1 jM^R^Fr 


dT/dy 


(3.18) 

(3.19) 


Since the discrete flux current density vectors h M , h yM , and he are in general dis- 

continuous across the boundary of a solution element, the integrations in equations (3.6) 
- (3.9) are understood to be carried out on the interior of the conservation element but 
immediately adjacent to the boundary S. 

We are now ready to proceed with the evaluation of each of the integrals in equations 
(3.6) - (3.9). To conserve space, we will only go through the details for equation (3.6). For 
equations (3.7) - (3.9) we simply give the final results. 

Since each S(CE(i,j )) is a simple closed curve in E 2 , the surface integration required 

in equation (3.6) can be converted into a line integration form [1, p. 14] . With ds = da n, 
where n is the outward unit normal to S(CE(i,j)) and da is the length of a surface element 
in E 2 , we have 

ds d = dy i — dx j, (3.20) 



and 


h M • ds = — pv dx + pu dy = g M ■ dr 


(3.21) 


where 

de f ( \ 

9m = (~pv,pu) 


and 

dr d = dx i + dy j. 


(3.22) 


(3.23) 


The line integration is taken to be positive in the counterclockwise sense. If we denote 
the vertices of an arbritrary conservation element CE(i y j ) by P, Q, R, and S as shown in 
Figure 2, we have 

h M ■ ds = <b 9 m - dr 

S(CE(i,j )) JPQRSij ~ 

= f [j(PQ)m + J(QR)m + J(RS)m + J(SP)m}.. (3.24) 


where [J(PQ) M ]i,j denotes the flux of h M through the line segment PQij , and similarly 


for J(QR) m , J(RS)m , and J(SP) M . We then have (omitting i,j subscripts) 


rQ 

J(PQ) m d — J - py dx + pu dy 
= / — pv dx + pu dy 

Jxi + 


= PV. 

hi -V ~ 


^ 2 
dx 


with y = y } + 


Ay 


(3.25) 


J(QR) m d = f - pv dx + pu 

Jq 

= / — py dx + pu dy 

hi + %*■ 


dy 


- - j 




2 


pu dy 


with x = x, — 


A x 


(3.26) 
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rS 

J(RS) m = / - pv dx + pudy 

Jr 

r*‘ + ±r 

= I — pv dx + pu dy 

/•*«+V 

— I pv dx with y = y } 

Jx 


Ay 

2 


where 
p u 


p v 


ASP ) 



tie/ 


f P 

I — pv dx + pu dy 

Js 


/■w + V 

I — pv dx + pu dy 


. A_iL 


pu dy 


with x = x* + 


A x 

~jT 


/^o, 0^0,0 “ 1 “ (/h, 0^0,0 H~Po,o^ i,o )(^ *^i) “I” (Po,i u o,o “h Po, 0 ^0,1) (?/ 

~\~ ( Po, 0 ^2,0 “i" P2,0 ^0,0 "f“ | 01 ,O «i,o) *^i) 

“I - {Po, o^i.i “ 1 “ Pi f i ^0,0 "I” Po, 1 U 1,0 “ 1 “ p\,0 U 0,1 ) ,a 'i)(y 3 /j) 

(Po,2 ^0,0 “h Po,0 ^0,2 “i” Po, 1 ii 0 , 1 ) (y Vj) “I” H. O. T. 

Po,0 ^0,0 “H (/^1,0 ^0,0 “ 1 “ ^0,0 ^1,0) *^t) (Po.l ^0,0 Po,0 ^0,l) (y 

“H i.Po ,0 i^ 2, 0 “ 1 “ P‘ 2 , 0 ^ 0,0 “h /* 1,0 ^1,0) *^t) 

+ (/^o,o ^i,i H~ pi,i ^o,o "h po, 1 ^1,0 “h pi (0 Vo.i) — ar«)(y !//) 

+ (po,2 v Q<0 + po,o Vo, 2 + po, 1 v 0 ,x) (y - yj ) 2 + H. O. T. 


(3.27) 


(3.28) 

Vj) 

(3.29) 

Vi) 

(3.30) 
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All polynomial terms in equations (3.29) and (3.30) that are higher than second order 
have been represented through the abbreviation “H. O. T.”. These terms will be neglected 
throughout the remainder of the paper to be consistent with the second-order expansions 
(3.1) - (3.5) (i.e., h M is understood to be a second-order Taylor series expansion, and 

—* — ♦ — * 
similarly for h X M<> hvM, and h E )• 

Carrying out the line integrations in equations (3.25) - (3.28), we obtain 


J(PQ)m — ^2 (^0,0 ^2,0 4 /^.O ^0,0 4 />1,0 ^l,o) (3.31) 


A A v r 

4 Ax yp 0 0 ^0,2 4 Po,2 ^o,o 4 Po,\ 




Ay 3 

J(QP) M = 12 ^°’° w °' 2 ^°’ 2 U °>° P°> 1 W m) 


(3.32) 


’ A, x 2 A.x 

Ay | — (p 0t 0 U 2,0 4 p 2,0 W 0,0 4 p 1,0 U \,o) 2 (Pl.O U 0( 0 4 pQ t 0 W l,o) 4 ^> 0,0 ^ 0,0 


Ax 3 

J(RS)m — — (p0,0 ^2,0 4 p2,0 ^0,0 4 pl,0 ^l,o) 


12 


(3.33) 


Ax 


Ay 


Ay 


^ (Po,0 ^0,2 4 Po, 2 Vq,0 4 P 0,1 ^0,l) 2 (/**<>, 1 4 p0 t 0 ^0,l) 4 po t o ^0,0 


Ay 3 

J {SP)m — ^2 (p°>° W 0,2 4 Po t 2 U 0,0 4 Po t i W 0,l) 

r Ax 2 Ax 

4 Ay — ( po,o u 2t o 4 p 2,0 u 0> o 4 pi,o ^i,o) 4 ^ (pi,o ^0,0 4 po,o ^1,0) 4 Po t 


(3.34) 


0 W 0,0 


For clarity, all i,j subscripts have been omitted in equations (3.31) - (3.34). By virtue 
of equations (3.6) and (3.24) we require that 

J(PQ) m + J(QR) m + J(RS)» + J(SP) M = 0 (3.35) 

which simply says that the total flux of h M out of S(CE(i,j)) is equal to zero. Applying 
this condition to equations (3.31) - (3.34), we obtain the mass flux conservation constraint 

Pl ,0 ^0,0 4“ p 0,0 M 1,0 “1“ pQ.l V 0,0 4” Po,0 V 0> 1 = 0. (3.36) 
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Imposing this condition, we obtain from equations (3.31) - (3.34) the following expressions 
for the normalized mass flux across the boundaries of the solution element SE(i,j ): 


JjPQ)*, 

Ax 


Ai' 

~n 


(Po ,0 ^ 2,0 d~ Pl ,0 U 0,0 4" pl ,0 Wl.o) 


+ 


A y\ 

^ \P 0,0 V 0,2 


~\~ Po, 2^0,0 Po, 1 Uo.l) 



0,0 


Po.o ^ 1 ,0 ) 


(3.37) 


+ p 0,0 ^0,0 


J(QR)m 

Ay 


Ay 2 

12 


(Po,0 ^0,2 d" Po,2 U 0,0 d - ^0,1 ^0,l) 


(3.38) 


Ax ^ / vt 

d (pO.O U 2,0 d- P2.0 W 0 ,0 d- Pl ,0 Ki.o) — “ 2 ” (^ 1.0 u o,o + Po,o u i,o) 




JjRSU 

Ax 


Ax 2 

~12 


(po,o ^2,0 d - p2,0 V 0,0 d - p\,o ^ 1 . 0 ) 


(3.39) 


Ay* , Ay 

d“ ( po,0 V 0,2 d" Po,2 V 0,0 d - P 0,1 Wo.l) d" ■ n \P>.0 U 0,0 


JjSPU 

Ay 


12 


(Po,0 ^0,2 d" Po,2 U 0,0 d” P 0,1 


(3.40) 


Ar 


d — (po,o U 2,0 d- ^2,0 u o,o d- Pl.o u i,o) d- ~n~ (P \,0 Uo,o d- Po,0 




Note that with these expressions for J(PQ) M , J(QR) M , J(RS) M , and J(SP) M , 
condition (3.35) is satisfied exactly. 

The mathematical procedure required to obtain the discrete flux conservation equa- 
tions corresponding to equations (3.7) - (3.9) is exactly analogous to that given above for 
equation (3.6). 

Corresponding to equation (3.7), we obtain the x-momentum flux conservation con- 
straint, 


1 2 

/>o,o(«i,o«o,od«o,iV) d- Pl.o - -j^-[(2u 0 .2 d- Wi.i) d- j(4«2,o 
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and the following expressions for the normalized flux of h XM across the boundaries of 
SE(i,j): 


JjPQh 

Ax 


(3.42) 


Ax 2 

~ 12 ~ 


^o,2 Po,o ^0,0 ^0,1 (/^0,1 ^0,0 *f /^0,0 ^0,l) 


/ Ay 2 

^2,0 Po,0 ^0,0 “h ^1,0 (Pi ,0 ^0,0 “I" p0,0 ^1,0 )j “f“ 

y ^ ' , 0 Vo,oU 0 ,i ~ -^^—(2 tr 0>2 + - -^-(«o,i + u,,o) + u 0 ,0 


+ 


P 0, 


Ax 


jXgg)j 

Ay 


Ay 


2 r 


+ 


12 
Ax 2 r 


^ 0,2 Pq ,0 ^ 0,0 4" ^ 0,1 (p 0 ,l ^ 0,0 4” Po.O ^ 0,1 ) 4” Po ,2 


Ax r 1 /n ,1 

4 — — |p 0 ,o w 0>0 u 01 — -^—(2 1/ 0>2 + w 1>x )J + 


^2,0 Po,0 ^0,0 4 - 0 (p 10 U 0 o -(- Po,0 ^1,0) 4 - Pi, c 

2 


(3.43) 


P °.0 3 Re ^ Ul, ° v oa) 


u„ 


J(QRh 

Ay 


J(RS) } 


Ax 


Ax 2 r 


12 


^2,0 Po,0 V 0,0 4* U l,0 (^*1,0 ^0,0 4" p0,0 ^1 ,o )j 4- 


Ay 2 


Ay 

2 


P 0, 


(3.44) 

^0,2 Po,0 V<>,0 4" ^0, 1 (Po,l ^0,0 4" Po,0 Vo, l ) 

J(RS) M 


^0,0 ^0,1 J-J (2 U 0 2 4" U,.l) («0,1 4- Ui o) ^0,0 

afij; J Rc-t 


Ax 


12 



ASP) 

Ay 


X M 


(3.45) 


Ay 2 


u 


Ax r 

~~2 


+ 


P 0,0 Wo,o Wo, i 


12 

Ax 2 r 


0,lp0,0 U 0,0 + W 0 ,l (Po.l u o,o + po.O Wo.l) + Po,2 


U20 Po,0 Wo,0 4 " W 2 , o (P 1,0 Woo 4 ” Po,0 W,.o) 4 " 7 ^ 2,0 

1 2 

(2 u 0 , 2 4- U| (1 )| 4" Po,o (2wi_ 0 W 0 ,,) + w 0 ,o 


Re L 


Jjspfl 

Ay 


The normalized mass flux expressions (3.37) - (3.40) have been used to simplify equa- 
tions (3.42) - (3.45). 

Corresponding to equation (3.8), we obtain the y-momentum flux conservation con- 
strain t, 


P 0, 


,o (wj,o w 0 , 0 4~ w 0 ,j w 0 ,o) 4- Po.l — j^-[(2w 2 ,o 4- U,,,) + ^(4 r„, 2 w,,,)] _ 0 (3.46) 


and the following expressions for the normalized flux of h Y \t across the boundaries of 

SE(i,j): 


(3.47) 


Ax 2 


Ay 

2 


+ 


p o,o Wq o W, o 


12 

Ay 2 

4 


1 


Re, 


J{PQ)ym 

Ax 

2,0 Po,0 Wo,0 4 " W, o (p,,0 W 0| 0 4 " Po.O Wi,o) + P'2,0 
Wo, 2 Po , o W 0 , o 4 - W 0 il (po,iW 0 , 0 +p 0 l 0 W 0 > i) 4 - Po, 

1 2 tn ^ J(FQ)i 

(2 u 2 ,o 4- Wj i) + p 0 ,o — (2w 0 ,i w 10 ) + w 0 ,0 ^ 


JjQRfl 

Ay 


(3.48) 


Ay 2 

12 


r 1 

Ax 2 r , x' 

<3 

o 

to 

o 

o 

c 

o 

o 

+ 

o 

o 

o 

£ 

o 

+ 

o 

o 

o 

+ ^ |^2,0 Po,0 ^0,0 “f" ^1,0 (^ 0,0 Ul,0+Pl,0 u o t o) 


Ax r 

~2 


Po,o W 0 ,o W,,o - ^-(2 w 2>0 4- «i,i)j — ^(“o,i 4- w, °) V 0 ,0 


JjQRfl 

Ay 
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J(RS) 

Ax 


YM 


(3.49) 


Ax 2 


+ 


12 

Ay 2 


^2,0 Po,0 Vq,0 d" ^1,0 (Pl,0 ^0,0 “I” Po,o ^i,o) d~ Pi, o 


+ 


Ay 


^0,2 Po ,0 Vo,0 "I" ^0,1 ( Po,l ^0,0 "I" p 0,0 Uo.l) "f" Po,2 

2 


1 1 2 

Pa , o I^o.o ^i,o (2 ®j,o "f Ui,i) 4" Po , o — n (2 Uo , i M i o ) 

AC^ J UilCr 


J(PS) 


M 


Ax 


j(gP)i 

Ay 


(3.50) 


Ay 2 

12 


i r 

W 0 ,2 Po ,0 ^0,0 “i“ ^0,1 (Po, 0 Wo, 1 H - " Po, 1 W 0 o)j ^ |^2, 0 Po, 0 Wq, o + ^1,0 (Po,0 Wi 0 "t"pl ( 0 W 0) o ) 

— \« uv _ 1 (2v -1 -u )] _ 1 ( u + v \ _i_ v J ( SP )m 

0,0 W 0i 0 ^ 1,0 \^ V 2,0 I ^ 1,1 / J V W 0,1 ■ ”l,o) I W 0 0 


Ay 


Corresponding to equation (3.9), we obtain the energy flux conservation constraint , 


Po.o («I.O + Uo,l) + 


Po,a (T,.o Uo o -J- T 0 l Uq o) 


Re, 


^7(7 - 1) 

^1,0 (^i,o d - ^0,1) d" ^o,i (^i,o *t ^0,1 ) d - 3 1,0 (2 U| 0 ^0,1) d - ^0,1 (2^o,i ^1,0)] 

2 


( 7 -l)M^e L Pr 


(^2,0 d- ^0,2) = 0 


(3.51) 


and the following expressions for the normalized flux of h E across the boundaries of 

SE(z,j): 
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J(PQ), 

Ax 


(3.52) 


Ax 5 

~12 

4 


| /Vo V ,0 ( 1,0 — + U 0,0 U 2 ,0 4 I’o.O V,o) + (p 1,0 V.O 4 Po,0 V.o) («1.0 ^ 0.0 4- Wl.O V,o) 

1 

[T 10 p Qt0 V l 0 + ^o,o(^Vo Po,o -^*1,0 P\ ,o)] “I" {v 2 0 p 0 0 V x q Pt o V 0 0 p 7 0 ) 


MZo 7(7 - 1 ) 


3 Re: 


•ko( 2 «o,,-Ui,o) 4 2 u,,o (u l ,i-« 2 ,o)] - ^-Ko(v,,o+tio,i) + «i.o( 2 u a ,o + «i.i)]} 


2 2 | 2 

+ — ^ Po,0 V 0,0 ( 1 2 1 I _ 1 £ 0,0 11 0,2"I" V 0.oV,2) 4 ( Po,l v o,o “I" Po,o V 0, 1 ) { U 0. 0 U 0 , 1 4 U 0,0 V 0,l ) 

[T 0i 2 po,0 V 0,0 4 T 0 'i (p 0 ,l V,0 4 Po,0 U 0.l)] 4 " (^0,2 Po,0 4 V 0,\ Po,l 4 V 0,0 P0.2) 


4 


1 


m £> 7(7 - 1 ) 

1 


Re L 


[u 0 .a(Vi,o+U 0 .i) + u o.> (2tt 0 ,a + »i.i)] “ K 2 (2 U 0> i -«..o ) 4 V 0l i (4 V 0 , 2 } 


— ^ /V 0 ^ 0, 0 ( V, 0 ^ 0 , 0 4 V|, 0 V, 0 ) 4 o Po,0 4 Pl,0 ^ 0 , 0 ) 4 y ^ ,0 Po,0 ^0,0 


3i?e, 


Re, 


^1,0 (2^1,0 V,l) 4 «o,o( 4«2,0 V,l) 

2 


V,o (uj.o 4 ^0,1) 4 V,o (2 ^2,0 4 ^1,1) 


( 7 - 1 )Ml 0 Re L P 




Re L 


^ 0,0 (^i,o 4 ^ 0,1 ) 4 3 v 0 ,o (2 v g 1 Ui 0 )j 4 v 0 ,o Po.o 


1 


(7 - 1 )Ml 0 Rt L Pr 


T 0>1 




JjPQh 

M^7(7-1)‘ U '^ Ax 


°) 
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mxh 

Ay 


.. , u l. i+ v o,i , .. , „ .. N , . .. w . .. , 

22 \ 2 ' ^0,2 I ^0,0 ^0,2 ) + (^0,1 ^0,0 + p 0,0 tio.l J (^0,1 ^0,0 + ^0,1 ^0,o) 

+ < y^ / y i) ^ 0, ° ^o,o(-^o,2 Po,o “I” -^0,1 Po p i)] + (^o p 2 Po.o + ^o,i Po p i + ^o p o ^0,2) 

2 2 
3^Re ^ 0,2 ^^ l *° ^°> 1 ) 2t/ 01 (ti lpl ^ 0 , 2 )] [ v o,2 (^i P o ~ I - ^ o, 1 ) “I - ^o.i(2u 0i2 +i; 1(1 )]} 


j^o.o U 0,0 ( 


u lo + v lo 


+ ti 0 o ^2,0 + ^o,o ^2 , 0 ) “ 1 “ (Pi,o ^o,o “I” Po, 0 Ui (0 ) (^0,0 ^1,0 + ^0,0 ^1,0) 


+ \/r*i { 1 \ PVo po,o u o,o + Ti.o (Pi,o ^o,o + Po,o ^i p o)] ” 1 " (^2,0 Po,o + ^i p o Pi p o “I - ^o,o £*2,0) 

M ool\l - l ) 

2 1 >1 

- g^-[u 2i0 ( 2 « 1 | 0 ~V 0 p i) + Ui ( 0 ( 4 , 2i0 ^ m )] - ^-ko(Vi,o + Uo,i) + ^1,0(2^2,0+141,1)] | 


^/?0,0 ^0,0 (^1,0 ^0,0 + t>i,o ^o p o) + (^1,0 Po,o + Pi,o^o p o) + 


- 1 ) 


T l|0 P 0,0 ^0,0 


^1,0 (2 U 10 ^ 0 , 1 ) + ^0,o(4t42,0 ^1,1 

°L *■ 


Re ^ l, ° ^ l, ° Wo ’ 1 ^ ^ 0, ° w m)J 


(7-1 )Ml Q Re L Pr 2 ’ 


jj ^0,0 (^1 ,o + ^o, 1 ) + 3^0,0 ( 2 l 4 i 0 ^o,i) + ^o p oPo p o 

XtCr L J 


(7 — 1 )Ml 0 Re L Pr 


( U lo+ V lo 1 rp \J(QR)m 

V 2 + M^ 7 ( 7 -l) io 'V Ay 
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J(RS)e 

Ax 


(3.54) 


Ax 2 

HT 

+ 


2 2 

|/Vo ^0,0 ( 1>0 -~ + Wo,o u 2,0 + ^0,0 i; 2 p o) + (pi,o v o,o + Po,o Vl.o) (Ui.o Mo.o + Vi.o Uo.o) 

[T 1>0 Po,o W,, 0 + W 0 , 0 (T 2 ,o Po,o 4* T,,o Pi.o)] + (w 2 .o Po,o + W, 0 p 10 4* W 0 0 p 2 ,o) 


M^ 7 (7 - 1) 


3i?e, 


■ w 2 , o (2 ' 


0,t — Ml.o) + 2i», i0 (Vi.i — «2,o)] — ^“[ U 2,0 (Vi,o+«0,l) + U li o(2t> J> o+«, i |)]| 


+ 


Ay 2 


^Po,0 W 0 ,o ( 


«0, l 4* v o,i 


+ 


4 L r —V 2 

1 


+ U 0,0 u 0,2 + W 0 , 0 Wo, 2) + (Po.l W 0 ,0 + Po,0 W 0 ,l ) (« 0,0 U 0,l + W 0 ,0 Wo,,) 


Ml 7(7 - 1) 


[T 0 , 2 Po,o W 0 ,d + To,, (p 0 ,l Wo.o + Po,o W 0 ,i)] + (w 0 ,2Po l0 + W 0 , 1 Po.l + W 0 ,0 P0.2) 


— — [«o,2 (w,,o + «o,i) + u o.i (2«o,a + u i,i)) ~ 
iic L 


3 Re, 


'[w 0 ,2 (2w 0 ,, u,,o) + W 0 ,, (4 u 0>J -«,.,)]} 


+ jft.O U 0,0 ( U l,0 U 0,0 7 ^1,0 Wo,o) 7" (WjoPo.o+Pi.oWo.o) 4" M 2 -y^-y \ ) T, 0 Po,0 W 0 ,o 


3 Re, 


u 


Re, 


1,0 (2 U, o W 0 ,i) 4" ^0,0 (4 ^2,0 W,,i) 

2 


Wi,o (w,, 0 + w 01 ) 4- w 0 , 0 (2 v 2 0 4- “1,1 )| ('y — 1)M 2 Re^Pr "^ 2, ° 1 


Re, 


u o,o (wi.o 4" w 0 ,i ) 4- 3 w 0 ,o (2 u 0 , , u i ,o)j 4- w 0 0 Po,o ^ — \)M 2 Re i^Pr^ 0 ' 


_ ^0,0 


Wo o + W 0 2 0 


4- 


J(RS), 
Ml 7(7-1)'“’^ Ax 


•°) 
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JjSPh 

Ay 


(3.55) 


Ay* 

12 


r ^i 2 H" t? 2 

Po,0 ^ 0,0 ( 2 ^ ^ 0,0 ^ 0,2 H“ ^ 0,0 ^ 0 , 2 ) "I" (p 0 , 1 ^ 0,0 ” 1 " Po,0 u 0,l)(u 0tl U 0 Q ~}~ i V 0>0 ) 


4 ” j ^2 / y(^'y 1 ^ 4 “ ^o,o(-^o,2 Po,o 4 ~ -^0,1 Po,i)] 4 ” (^o,2Po,o 4 “ ^0,1 Po.i 4 ~ ^0,0 Po,2) 

2 y 

~ 2 ^”[ w o ,2 (2w 1 | 0 -v 0 l ) + 2 u 01 (u lfl — 1 > 0(2 )] — TJg [^o,2 (^i,o “h w 01 ) + Vo,i( 2 iio #a +t; lll )]} 


+ 


Ax 2 r 


^Po.O ^0,0 ( 


y? 1 ^2 
U l,0 ' U l,0 


4 " ^ 0,0 ^2,0 “ 1 “ ^ 0,0 ^ 2 ,o) “I” (Pl ,0 ^ 0,0 4 " p 0,0 ^l,o) (^ 0,0 W t)0 “ 1 “ ^ 0> o ^1,0) 


4“ M 2 7(7 — 1) ^ 2>0 ^ 0, ° Wt>, ° ■^ 1)0 (^ l -° W °-° 4” Po ,0 Wl *°)] + ( w 2,0 po.o "f ^ 1,0 Pl ,0 + W 0( o P 2 ,o) 

2 1 ^ 
' [ w 2,0 (2 ZX 1)0 — V 0(1 ) + w 1(0 (4u 2 i o~Ui,i)] — [ v 2,0 (^ 1,0 ~l~ W 0,l ) + t^io (2 t?2,0 4" Wl,l )] J 


3-Re, 


Ax r 1 

4 ” r> 1 Po.o w o,o («i,o w o,o 4 “ ^1,0 ^o,o) 4 “ (^i.opo.o 4 “ Pi,o ^o,o) 4 - ~rjo 7 TV Ti.o Po,0 0 
2 l M 2 )7 ( 7 -l) 

^1,0 (2u 1(0 ^O p 1 ) 4” ^0,0 ^2,0 ^l,l)j 


3Re, 


Re, 


^ 1,0 (^ 1,0 4- W 0 ,l ) 4“ ^0,o(2'^2 ) o4~W 11 )j 


( 7 - l)M^Re L Rr 2 


°} 


Jfc J»o,o(»i,. + tt, lI ) + ? u o,o( 2 tx li0 t; 0il )] + « 0 i 0 p 0i0 T lf0 


, ( ujo+vlo 1 T \ J(5P) 

V 2 + M^ 7 ( 7 -l) io ’V Ay 


M 


— r ^ ^ 

With these expressions for the flux of h XM , h YM , and h E across a cell boundary, the 
following conditions are satisfied on every solution element. 


J(PQ)xm + J(QR) xm + J(RS) xm + J(SP) XM = 0 


(3.56) 


J(PQ)ym ~b J(QR) ym + J(RS) ym + J(SP) ym = 0 


(3.57) 


J(PQ) e + J{QR) e + J(RS) b + J(SP) E = 0 


(3.58) 
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Before concluding this section, one final comment is required. For compressible flows, 
it will be necessary to assume an equation of state in order to close the system of equations 
(2.11) - (2.14) and their discrete counterpart (3.6) - (3.9). Since we have assumed that the 
fluid is a perfect gas, we may take (using dimensional quantities) 

p* = p* RT* (3.59) 


and its discrete analogue 


p * = p*RT*. 


(3.60) 


In dimensionless form, we then have 


P = 


7 Ml 


PT = 


7*C 


P 0 


,oTo,o + (pl.o ^0,0 + /Vo ^I.o) ( x — x i) (po,i Tq.o + Po.o ^0,1 ) (j/ Vj) 


“h (/?{), 0 ^2,0 “1" P'2,0 T 0>0 pl,0 ri i0 ) ^») 


(3.61) 


d (po.o T\,\ d Pi, i To.o d Po,\ T 1>0 d pi,o To , i ) (® •*•>)(?/ Vj ) 


so that 


d (Po,2 T 0 , 0 -(- Po,0 F 0 2 d Po.l T 0 ,i) (t/ 2/jf ) | 


rfe/ 1 

Po '° - 7 M*/ 0l ° 0l ° 


(3.62) 


Po.l — „ ,9 (Po.,T„,o+Po,oT 0 ,,) 


(3.63) 


Pl,0 — , r9 (Pl,0 T 0 , 0 + Po.O T, 0 ) 

7^5o 


(3.64) 


Po 2 d — —T7T (Po,2 T 0 ,0 d Po,oT 0 ' 2 + Po,l ^O.l) 


(3.65) 


Pi,i — — . , 2 (Po.o Ti.i d Pi,i To.o d Po.i Ti.o d Pi,o T 0>1 ) 

iMoo 


(3.66) 


P 2 , 0 — s ,9 (^0.0 ^2.0 d P2,0 T 0 ,0 + Pl,0 T,,o). 

7^5» 


(3.67) 
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IV. Application to Incompressible Channel Flows 


Our formulation thus fax has been for the compressible Navier-Stokes equations. Ap- 
plication to the incompressible equations can be done very simply as a special case. 

When the density is constant, the gradient terms p 0 ,n Pi.o, Po.i, Pi.i, and p 2 0 in the 
Taylor series expansion (3.1) are all zero, and the constant term p 0i0 may be set equal to 
one. Equations (3.31) - (3.58) then apply to incompressible flow, and equations (3.59) - 
(3.67) are no longer applicable. 

Consider the channel geometry and mesh shown in Figure 3. On each solution element 
SE(i,j ), there are 4 unknown discrete variables, u,v,p , and T, each with 6 unknown 

coefficients, for a total of 24 unknowns per cell. However, the mixed coefficient terms p, , 
and T,.! do not appear in any of the discrete equations (3.37) - (3.55). We will assume 
that p, , = T, 4 = 0 as a result. In addition, the energy equation decouples from the 
continuity and momentum equations, so that the 17 unknowns for the two components of 
velocity and pressure may be solved for independently of the five remaining unknowns for 
the temperature. 

Let Nj denote the number of solution elements from the lower wall to the upper wall, 
and let N{ denote the number of solution elements in the downstream direction. Then 
the total number of unknowns for the velocity and pressure is 17 N,Nj. We thus require 
17 NiNj conditions to close the system. 

The flux conservation constraints (3.36), (3.41), and (3.46), which ensure local flux 
conservation on each cell, provide the starting point of the present formulation. For in- 
compressible flow, these become 


u i,o "b ^o,» 


= 0 




(4.1) 


^i,o ^o,o “l - ^o,o “1“ Pi.o [(2tXo 2 "h ^i,i) “I - g (4u 2>0 Wi,i)]j — 0 (4.2) 


1 2 

v 1}O u 0t0 + v 0jl v 0t0 + p 0(l - “ — [(2 v 2>0 + u 1(l ) + - (4 v 0>2 - u ltl )\ 

Ib€- £ tJ 


. . = o (4.3) 


for j = 1 , . . . , N j and i = 1 , . .. , N t , for a total of 3 NiNj conditions. 


To ensure global flux conservation, mass and momentum fluxes must be balanced at 
cell interfaces. Consider the vertical interface between SE(i,j) and SE(i + l,j), as shown 
in Figure 3. Using the cell orientation PQRS previously introduced, the mass flux leaving 
SE(i,j ) through the vertical interface is given by [J(5P) M ]ij. On the other hand, the 
mass flux leaving SE(i + 1, j) through this same interface is given by [J(Qil) M ]i+i ,j (and 
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similarly for the momentum fluxes). Conservation of mass (and momentum) requires that 
the sum of the two be zero. Balancing fluxes across each vertical interface in the mesh, we 
thus require that 


' j(SP) u ] 

Ay 

L J »,j 

~ J(5P)xJ 

Ay 

L J *,j 

J(SP) Y M 1 [ AQR) ym 1 __ o (4 6) 

Ay . . + Ay " 

L J ij L J i+1,j 

for j = 1 , Nj and i = 1, ..., TV; - 1, for a total of 3 Nj(Ni - 1) = 37V, TV, - 3 Nj conditions. 
Similarly, the fluxes through the horizontal interface between SE(i,j ) and SE(i,j + 1) 
must also balance. We thus require that 

'jjPQU 

Ax 

'J(PQ)XM 

Ax 

' J(PQ)y M 

Ax 

for j = 1, Nj - 1 and i = 1, ..., JV,-, for a total of 3N t (Nj - 1) = 3 N { Nj - 3 N, conditions. 
Equations (4.1) - (4.9) thus represent a total of 9 N,Nj - 3 Ni - 3 Nj conditions. 

In addition to the above requirements, there are discrete boundary conditions that 
must be satisfied. For each solution element along a wall boundary, we must require as 
a minimum that there be no mass flux through the wall. This can be accomplished by 
setting the integrated mass flux at the wall to zero, or by setting the v velocity component 
identically to zero on the wall boundary. The second approach requires three conditions, 
whereas the first approach requires only one. In this paper, we use the first approach. 
Note that the zero mass flux condition takes the place of a wall boundary condition for 
the v velocity component. For the u velocity, we simply set u = 0 at the midpoint of the 

wall face of each cell. This leads to the following 47V, conditions: 
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0 


(4.10) 


rAy‘ ay 1 

^ Vo, 2 2 W ®.» "b ^O.OJ . J — 


Ay 

Ay 2 Ay 

r 4 v 0 , 2 

' Ay 2 , Ay 
. 4 


'Ax 2 j, 

^2 V2 '° "b ^ v o,2 2 w °- 1 ~b *V®J . j = 0 


, Ay 1 

Vo, 2 T „ V 0 j T v 0 0 — 

2 J i,Nj 


(4.11) 

(4.12) 



'Ax 2 

12 

Ay 2 Ay ] 

V 2 , 0 + . Vo, 2 + - V 0ll + V 0 o I = 0 

4 2 Ji \,Nj 

(4.13) 

Boundary conditions must also be specified at the inlet and exit. 

At the upstream 

boundary, we specify plug flow inlet conditions: 




■ At^ /\ T 1 

4 U2 ’° 2~ Ul ’° + u ° 0 ] lj - = 

(4.14) 



Ax 2 Ax ] 

. V 2 o n V,,o + V(, o = 0 

4 2 j 1,7 

(4.15) 


Downstream, we specify the pressure: 


P 0,0 


N it j 


= Pe- 


(4.16) 


The total number of conditions is then increased to 9 N t Nj -f N t . 

Equations (4.1) - (4.16) are the minimal conditions that must be satisfied and are the 
starting point of our numerical formulation. Note that conditions (4.1) - (4.9) are integral 
relationships, and that these conditions ensure that fluxes of mass and momentum are 
conserved on every cell and every union of cells in the mesh. Having ensured that these 
fundamental requirements are satisfied, we may now turn to the differential form of the 
discrete equations to introduce additional conditions to close the system. We may thus 
require that 

V • (u, v) = V • h M =0 (4.17) 


V 

( V 2 + p - T xx , UV - 

- I*v) = V- 

’ fyxM 

= 0 

(4.18) 

— # 
V 

•(uv - T X y, V 2 + P - 

II 

> 

t 

hy Af 

= 0 

(4.19) 


on each solution element. (These relationships are already satisfied at the cell center by 
virtue of the flux conservation constraints (4.1) - ( 4.3), but only to first order throughout 
the rest of the solution element.) This leads to the following 6 N{Nj conditions: 
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2 u 20 + v ltl \ = 0 


2t) 0| 2 “I - ^1,1 ^ 

J m 


(4.21) 


2(2«o,„« J ,o + Wi,o + Pi, o) + Ui,i v o,o + VijUoo + U,,«Uo,i + u o,i w l,°].^ — 0 (4.22) 


l 0i 0 Uo.J + V 0 ,0 U 0, 2 + V 0 l U o l + Ul,l«0.0 + «1,0«0,I . . 

J I , J 


(4.23) 


2 (2 + vl, + p,.,) + »ot. l0 + = 0 < 4 ' 24 ) 


u 00 v 20 + u 0 0 u 2 0 + u i,o u i,o 4" v i,i v o,o 4- v 10 v 0- i 0 (4. 2d) 

j i,j 

The total number of conditions is then increased to l5N{Nj + N, . If we assume that 
p 2 o = 0 in view of the fact that channel flow is dominated by gradients in the cross-flow 
direction, then the total number of conditions required is reduced to lQN t Nj. We thus 
need N 2 Nj - N, — Ni(Nj - 1) conditions to close the system. Since Ni{N } - 1) is 
the number of horizontal cell interfaces in the mesh, this suggests imposing a condition at 
each horizontal cell interface. An obvious choice would be to require that the u velocity 
component be continuous in the cross-stream direction. We then have as our final condition 


Ay 2 Ay l rAy 2 Ay 

— Wo, 2 — ~7T U M 4" «0,0 . . — — J~ U 0,2 4- 9 «0,1 4- U 0,0 


(4.2C) 


for j = 1 - 1 and i = 1, ...» iVf. 

We now summarize the present implicit scheme. Using (4.1), (4.20), and (4.21) to 
eliminate v 0>l , u, ,, and tz M , and dropping the p 2 , 0 term, there are 13 unknowns per cell, 
and the Taylor series expansions become 


u(x,y,i,j ) d = u 20 (x — Xi) 2 — 2v 02 (x — Xi)(y yj) + u 02 (y yj) 
4“ tii,o(® ®«) 4- u 0 ,i(y yj) 4- «o,o 


(4.27) 


v(x,y\i,j ) = v 2 ' 0 (x — Xi) 2 — 2 u 2f0 (x — £j)(y y>) 4- u 0 2 (y yj) 

4* v l0 (x — ^i) — «i,o(y — yj) 4- v o,o 


(4.28) 


p(x,y,i,j ) = Po, 2 (y - V}) 2 + Pi,o(z — £i) 4-p 0 ,i(y y;) +Po,o (4.29) 
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Equations (4.1) - (4.3), which are the flux conservation constraints corresponding to 
equations (3.6) - (3.8), become 


U l ,0 + 


t>o,i = 
J •>> 


(4.30) 


*^1,0 ^0,0 4" ^0,1 ^0,0 "b Pl,0 r, (^J.O ~b ^ 0 , 2 )] — 0 

L ne L J i,j 

^1,0 ^0,0 ^1,0 ^0,0 + Po,l (^2,0 "b ^ 0 , 2 )] = 0 

L ttCi J i,j 

Equations (4.4) - (4.6), which are streamwise interface conditions, are given by 

J(SP)m 


(4.31) 

(4.32) 


Ay 


+ 


Ay 


J t,J 

Ax 2 


J(QR)» 

Ay 
Ax 


*+l J 


10 ^0,2 4“ A ^2,0 4" 0 U U o 4“ ^0,0 

12 4 2 




Ay 


Ax 2 


Ax 


10 ^0,2 "b , U 2,0 0 U l,0 "b u 0,0 

12 4 2 


= 0 




(4.33) 


4(5P) xm ‘ 

+ 

Ay 



J(QR), 

Ay 


*+i,j 


Ay 2 Ax 2 

(2 u 0 2 u 00 + ^0,1 "b £*0,2) "b ~ ( 2 u 2 0 u 0 0 -b rt, 0 ) 


12 


Ax , 


(4.34) 


[t?o,0^0,l ^0,0 ^1,0 


2 2 

( W 0,2 ““ W 2,o)] 4 " Po,0 — Wl,0 4 " U 0 0 


* J 


/\ 7 * ^ 

(2u 0> 2^o,o 4“ xx 01 4“ £*0,2) 4~ t (2t/2, 0^0,0 4~ w 10 ) 


12 


Ax r 2 N1 2 

4- 1 ^ 0,0 ^0,1 U 0,0 W l,0 (^0,2 ^2,0/] 4“ Po,0 T-) 

2i iiCi 116^ 


0 


u 


i.o 4 - u 0Q 


•+i,j 
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JjSPh 

Ay 


4- 


M 


12 


J(QR) 

Ay 


Ax 2 


YM 


(4.35) 




(^0,2 u o,o + W 0, 2^0,0 “ W 10 U 01 ) + — \V 2i oU 0)0 + W 2, 0^0,0 + V l(0 Ui.oJ 


Ax 2 1 
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Equations (4.7) - (4.9) and (4.26), which are cross-stream interface conditions, become 
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The boundary conditions at the lower and upper walls are 
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Conditions (4.20) - (4.25), which are derived from the differential form of the governing 
equations, are given by 


^i i + 2 u 2(0 


u xl + 2 v 0t 




hj 


= 0 


= 0 


2 (^ 0,0 ^ 2,0 ^ 0, 0 ^ 0,2 ) “H ^ 1,0 * 4 " Wo,l ^ 1.0 


= 0 




^0,0 ^0,2 ^ 


0,0 ^0,2 — 
J *ii 


2 (V 0t0 V 0 ,2 U 0,0 W 2,0 4” ^ 0 , 2 ) 4” ^i ( 0 4” U 0,l U I,0 


- 0 




U 0, 0 ^ 2, 0 ^ 0, 0 ^ 2, 0 


= 0 




and the upstream and downstream boundary conditions are 
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(4.52) 


Equations (4.30) - (4.52) are a coupled system of second-order polynomial equations 
that take the form 


F(X) = 0. (4.53) 

The unknowns p 10 , p 0> , and p 0 2 may be eliminated using (4.31), (4.32) and (4.48), respec- 
tively, so that there are a net 10 unknowns per cell that must be explicity solved for. The 
total number of equations in the nonlinear system is 10 N{ Nj. 

The Jacobian matrix associated with equation (4.53) is block tridiagonal, with block 
sizes equal to 10 Nj. If the equations in (4.53) are arranged appropriately, the Jacobian 
matrix has the structure shown in Figure 4. Assuming an initial solution iterate .V°, 
Newton’s method for this system takes the following form: 
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dF-\" 

M. 


(4.54) 


A n X = — F(X n ) 
where 


x n+1 = X n + A n X (4.55) 

The Jacobian matrix shown in Figure 4 has a highly diagonalized block structure. 
Ninety per cent of the lower diagonal blocks, and 80% of the upper diagonal blocks, are 
zero. In addition, there are no equations that span all three blocks in a block row of the 
matrix. If pivoting is employed during the elimination process, the upper diagonal blocks 
will fill in, and the storage required is j~N{N 2 — “TV 2 where N = ION, is the block 
size. However, if the equations that are located entirely within the diagonal block are 
pivoted separately from the equations that span two blocks, there will be no fill-in, and 
the storage required is only jf N{N 2 — yhjV 2 . I n the former case, the iteration matrix is 
nearly block bi-diagonal. In the latter case, the iteration matrix is nearly block diagonal. 
For both cases, the Jacobian matrix retains a highly diagonalized block structure during 
the elimination process by virtue of the fact that none of the discrete equations span more 
than two blocks. 

The computational work required is significantly reduced in either case over the situa- 
tion in which the elimination must be carried out completely in all of the blocks. Calcula- 
tions performed on a Cray YMP have shown that, for the nearly block bi-diagonal iteration 
matrix, the computational work required per iteration is reduced by about 47%. For the 
nearly block diagonal iteration matrix, calculations have shown that the computational 
work required is reduced by about 60%. 
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V. Numerical Results 


In this section numerical results are presented from calculations of developing chan 
nel flows with Reynolds numbers of 100, 1000, and 2000. Our concern here is not with 
quantitative validation, but rather with demonstrating the overall features of the present 
scheme. Validation will be considered in a future paper. 

For plug flow inlet conditions, the incoming velocity profile may be arbitrarily speci- 
fied. We consider for convenience velocity profiles of the form 

ui(y) = C (- - y 2 )" - (5-1) 

As n increases, the inlet velocity profile approaches a top-hat shape. In this paper we take 
n equal to 6. (See Figure 5.) Calculations with n values ranging from 4 to 10 have also 
been performed but will not be presented in this paper. The constant C is chosen so that 
the integrated mass flux at the inlet is equal to one. 

One of the major goals of the present work is to demonstrate the feasibility of solv- 
ing developing viscous flows on coarsely spaced uniform grids. To that end, numerous 
calculations have been performed of developing channel flows on grids using from 4 to 12 
cells per channel width. Our results indicate that as few as 6 cells across the channel 
may be sufficient to resolve the developing boundary layer, depending on the Reynolds 
number. For each Reynolds number considered, our approach was to repeatedly solve for 
the developing flow on successively finer grids until the predicted boundary layer thickness 
remained unchanged from the previous calculation. In the results that follow, we present 
results from the two finest grids that were used for each Reynolds number. We should 
point out that, due to symmetry, the number of independent cells to resolve the flow field 
is half the number of cells per channel width. 

The computational grids used in the present study are shown in Figures G - 9. In each 
case the exit boundary is 11 channel heights downstream. Figures 10 - 12, 13 - 15, and 16 - 
18 present results for Reynolds numbers of 100, 1000, and 2000, respectively. The predicted 
streamwise velocity profile is shown at 1, 3, 5, and 11 channel heights downstream. Each 
figure also shows the inlet velocity profile and the fully developed analytical solution. 

Figures 10 - 11 clearly indicate that at a Reynolds number of 100, the flow becomes 
fully developed in only a few channel heights downstream. At x = 5, the predicted stream- 
wise velocity along the centerline differs from the fully developed solution by .35 %. At 
x = 11, the predicted profile and the fully developed solution agree to a minimum of 3 
decimal places everywhere. Figure 12 combines the results from Figures 10 and 11. The 
results show clearly that 6 cells per channel width are as adequate as 8 cells to resolve the 
developing boundary layer. The CPU times required for these cases were 2.4 and 5.0 sec- 
onds, respectively, on a Cray YMP. The solutions were converged to a maximum residual 
error of 10 -4 , starting from an initial guess of uniform flow. 

At the higher Reynolds numbers of 1000 and 2000, the boundary layer is thinner 
and develops much more slowly. Consequently, more than 6 cells per channel width are 
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required to resolve the developing flow. The numerical results in Figures 13 - 18 indicate 
that 8 cells across the channel are sufficient to resolve the boundary layer at a Reynolds 
number of 1000, and 10 cells are sufficient at a Reynolds number of 2000. The CPU times 
ranged from 6.1 to 9.9 seconds for the Reynolds number 1000 results, and from 11.1 to 
16.6 seconds for the Reynolds number 2000 results. 

The numerical results in Figures 13 - 18 also show the occurence of non-physical 
oscillations in the velocity profile. These occur as a lingering effect of the singularity at 
the channel inlet. At the lower Reynolds number of 100, the naturally occuring viscocity is 
sufficient to quickly damp out the effects of the singularity. However, at the higher Reynolds 
numbers, this is no longer the case, and the oscillations persist further downstream. We 
should point out that a more physically realistic inlet velocity profile would reduce the 
strength of the singularity and the resulting non-physical oscillations. 


Conclusion 


In this paper we have presented a new flux conserving numerical scheme for solving 
the two-dimensional, steady Navier-Stokes equations. There are numerous advantages to 
the scheme we have developed. First, fluxes of mass, momentum, and energy are conserved 
on every cell and every union of cells in a computational mesh which has been used to 
discretize a flow field. Second, fluxes are balanced at cell interfaces without the use of 
interpolation, extrapolation, or flux limiters. Third, the discrete solution obtained on each 
cell is a functional solution of both the integral and differential form of the Navier-Stokes 
equations. Fourth, the present scheme is highly localized, and concentrates most of the 
information on a local cell. This results in a nearly block diagonal Jacobian matrix with 
minimal solution times and storage requirements. Fifth, as shown above, the present 
approach offers the potential to solve developing viscous flows on coarsely spaced, uniform 
grids. Finally, the approach we have developed provides a unified treatment of the discrete 
dependent variables and their derivatives. All are treated as unknowns together to be 
solved for through simulating local and global flux conservation — i.e., physics. 
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Appendix 


Analytical Solution for Fully Developed Flow 

The analytical solution for the fully developed velocity is obtained by solving equation 
(2.2). For flow in an infinitely long straight channel, this equation reduces to 


cPu dp 

= tL dx 


(A.l) 


where ^ is the nondimensional constant pressure gradient, and u must satisfy the bound- 
ary conditions 

ti(-i) = 0 (A. 2a) 

u(i) = 0. (A. 2b) 

The solution to equation (A.l) and boundary conditions (A. 2) is 


Re L dp 2 

u(y) = — T x (v 


(A.3) 


For unit nondimensional velocity at the channel inlet, conservation of mass requires that 

(A.4) 


dp _ J2_ 

dx Re i 


The analytical solution for the nondimensional fully developed velocity is thus given by 

«(y) = - 6 (y 2 - ?)• ( A - 5 ) 

When the flow is fully developed, equations (4.30) - (4.52) reduce to a linear set of 
equations which can be solved analytically. For the fully developed case, the discrete 
equations reduce to 
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for j = 1 - 1 


r Ay 2 

A ^0,2 

Ay 

Uq 1 + U 0 o 

= 0 

(A. 9) 
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Equation (A. 6), which is the x-momentum flux conservation constraint, immediately 
determines the value of u 0 2 on each solution element. We have 


U 0>2 

J *,} 


Re i Re i dp 

~r p '° = ~2 ~Tx 


Rcl 

~ Y ~ 



(A. 11) 


since p, 0 is the known pressure gradient There are thus only two unknowns per solution 
element, u 01 and u 0 0 . If we consider the special case of only one solution element for the 
entire channel width so that Nj = 1, equations (A. 9) and (A. 10) form a system of two 
equations in the two unknowns u 0 ,i and u 0 0 . The solution of this system is given by 


Thus, 


«0,1 = 0, u 00 = §. 

(A. 12) 

u(y) = —6 (y 2 - i) = u(y) 

(A. 13) 


so that the analytical solution is recovered. 

In general, it can be shown that the solution to the system of discrete equations 
(A. 6) - (A. 10) is given by 
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(A. 15) 
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(A. 17) 


where 


y } = - i + O' - i)Ay. 


The analytical solution is thus recovered on each solution element. 


(A.i8) 
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A re-examination of equation (A. 6) indicates that 


Since 


we have 


2 ^0,2 — 1 
l i,j ax 

(A.19) 

d 2 u/dy 2 = 2 [u 0>2 

(A. 20) 

L j 


d 2 u/dy 2 = Re L ^ 
ax 

(A.21) 


which corresponds exactly to equation (A.l). The ^-momentum flux conservation con- 
straint for the discrete solution thus reproduces the governing differential equation for 
fully developed channel flow. 
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Figure 4. Jacobian Matrix Structure. 
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Figure 6. 39 X 6 Computational Grid. 
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Figure 10. Predicted streamwise velocity profile at 1, 3, 5, and 11 
channel heights downstream. Re = 100. 39 X 6 grid. 
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Figure 1 1 . Predicted streamwise velocity profile at 1 , 3, 5, and 1 1 
channel heights downstream. Re = 100. 39 X 8 grid. 
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Figure 12. Comparison of predicted streamwise velocity at 1, 3, and 11 channel heights 
downstream from calculations on 39 X 6 and 39 X 8 grids. Re = 100. 
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Figure 13. Predicted streamwise velocity profile at 1, 3, 5, and 11 
channel heights downstream. Re = 1000. 39 X 8 grid. 
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Figure 14. Predicted streamwise velocity profile at 1, 3, 5, and 11 

channel heights downstream. Re = 1000. 39 X 10 grid. 



Figure 15. Comparison of predicted streamwise velocity at 1, 5, and 11 channel heights 
downstream from calculations on 39 X 8 and 39 X 10 grids. Re = 1000. 
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Figure 16. Predicted streamwise velocity profile at 1, 3, 5, and 11 

channel heights downstream. Re = 2000. 39 X 10 grid. 
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Figure 17. Predicted streamwise velocity profile at 1, 3, 5, and 11 

channel heights downstream. Re = 2000. 39 X 1 2 grid 


Figure 18. Comparison of predicted streamwise velocity at 1, 5, and 11 channel heights 
downstream from calculations on 39 X 10 and 39 X 12 grids. Re = 2000. 
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